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1.  INTRODUCTION  AND  SUMMARY 


A numerical  procedure,  known  as  the  Time  Dependent  Computation  (TDC)  com- 
puter program  is  being  developed  which  determines  steady  state  solutions  of  the 
turbulent,  two-dimensional  Navier-Stokes  equations  as  the  asymptotic  large-time 
limit  of  the  full  unsteady  equations.  The  hyperbolic  character  of  the  unsteady 
equations  offers  the  advantage  of  solving  an  initial  and  boundary  value  problem 
which  is  well  suited  for  a digital  computer.  In  addition  to  providing  steady 
state  solutions,  accurate  transient  results  are  also  obtained  thus  adding  to 
the  program's  versatility. 

In  Section  2,  the  formulation  is  presented  for  time  dependent  calculations 
of  inviscid  and  laminar  flow  using  an  analytic  coordinate  transformation  which 
maps  the  space  between  the  body  and  the  initial  outgoing  wave  front  onto  a 
computational  plane.  A coordinate  transformation  is  developed  for  a symmetrical 
Joukowsky  airfoil  as  an  example.  Accurate  boundary  conditions  are  derived  for 
the  general  case;  a special  treatment  of  the  singularities  induced  by  the 
example  transformation  at  the  cusped  trailing  edge  is  also  presented. 

In  Section  3,  the  formulation  is  extended  to  time  dependent  calculation 
of  turbulent,  viscous  flow  by  incorporating  the  Extended  Energy-Dissipation  tur- 
bulence model  of  Spalding  into  the  equations  of  motion.  The  fine  mesh  spacings 
required  to  compute  flow  properties  across  the  laminar  sublayer  are  avoided  by 
developing  approximate  surface  boundary  conditions  based  on  the  "law-of-the- 
wall"  relations. 

Results  for  inviscid,  subcritical  calculations  demonstrating  convergence 
to  the  steady  state  solution  are  presented  in  Section  4.  Some  results  for 
completely  laminar  flow  and  laminar  flow  with  transition  to  turbulence  are 
presented  in  Section  5 for  the  case  of  a symmetric  Joukowsky  airfoil. 
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2.  TIME  DEPENDENT  COMPUTATION  (TDC)  PROGRAM  DESCRIPTION 

The  TDC  program  formulation  is  modeled  after  the  approach  first  pro- 
posed by  Moretti  (Reference  1),  with  particular  emphasis  placed  upon  proper 
physical  representation  of  the  flow  field.  The  flow  is  produced  by  accelerat- 
ing any  arbitrary  body  smoothly  from  rest  to  a constant  velocity  in  a fluid 
whose  motion  may  be  arbitrarily  prescribed.  The  numerical  computations  are 
performed  within  the  region  (physical  plane)  bounded  by  the  body  and  the 
surrounding  initial  perturbation  wave  which  propagates  outward  from  the  body 
at  the  start  of  its  motion.  This  allows  initial  conditions  consistent  with 
the  physical  situation  to  be  properly  imposed  along  with  the  enforcement  of 
correct  far-field  boundary  conditions  as  the  solution  develops  with  time. 

Correct  representation  of  both  initial  and  boundary  conditions  is  necessary 
to  ensure  convergence  to  the  correct  steady  state  solution. 

The  annular  shaped  physical  plane  surrounding  the  body  is  mapped  onto  a 
rectangular  computational  plane  by  means  of  a series  of  coordinate  mappings 
and  stretching  transformations.  This  is  illustrated  simply  in  Figure  1. 

These  transformations  are  chosen  such  that  uniform  grid  spacing  in  the  compu- 
tational plane  corresponds  to  a highly  non-uniform  grid  arrangement  in  the 
physical  plane.  This  provides  for  fine  coordinate  mesh  spacing  in  those 
regions  where  gradients  are  large,  such  as  boundary  layers,  wake  regions,  and 
separated  layers  and  coarse  mesh  spacing  where  gradients  are  small.  As  a 
result,  economic  computation  times  are  achieved  by  reducing  the  total  number 
of  points  required  to  adequately  describe  a given  flow  field. 

A steady  flow  solution  is  first  realized  in  a narrow  region  immediately 
surrounding  the  body  at  some  time  after  the  body  has  accelerated  to  uniform 
velocity.  The  size  of  this  region  continuously  grows  with  subsequent  computa- 
tions until  a steady  state  region  of  acceptable  size  is  attained,  at  which  point 
the  calculations  may  be  terminated. 

2 . 1 Navler-Stokes  Equations 

For  two-dimensional  problems,  the  Navier-Stokes  equations  are  most  efficiently 
solved  in  polar  coordinates  r and  6,  as  can  be  seen  by  examining  Figure  1.  Speci- 
fically, for  a perfect  gas  with  constant  viscosity  and  thermal  conductivity, 
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FIGURE  1.  COORDINATE  TRANSFORMATION 
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The  symbols  used  in  these  equations  are  defined  as  follows: 
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P 

Q 
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U 

V 
p 
r 
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<J> 

P 
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M 
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Re 

Pr 

P„ 
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- Reference  length 

- Log  (p/pj 

- Non-dimensional  temperature,  e 

- Log  (p/p  ) 


P-R 


Radial  velocity  non-dimensionalized  by  Vp/p 

r OO  00 


- Angular  velocity  non-dimensionalized  by  Vp J 

- Static  pressure 

- Radial  coordinate  non-dimensionalized  by  L 


- Time  non-dimensionalized  by  L\/o  /p 

- Dissipation  function 

- Mass  density 

- Polar  angle 

- Ratio  of  specific  heats 

- Free  stream  Mach  number 

- Free  stream  Reynolds  number 

- Prandtl  number 

- Free  stream  static  pressure 

- Free  stream  mass  density 


U - Longitudinal  free  stream  velocity  component 

V - Vertical  free  stream  velocity  component 


2.2  Transformation  Onto  the  Computational  Plane 


The  transformation  mapping  the  physical  plane  (r,6)  onto  the  computational 
plane  (X,Y)  has  the  form 
t -*•  T 

r - Y (t , r , 6)  (0  < Y < 1) 

9 -*  X(t,r,9) 

F(t,r,9)  - F (T, X, Y) 

where  F represents  any  dependent  variable,  P,  R,  U,  or  V.  Derivatives  of  F 
are  obtained  in  a straightforward  manner  from  the  above  transformation  func- 
tions. For  example, 

3t  3T  3X  3t  3Y  3t 
3F_3F3X3F3Y 

3r  3X  3r  3Y  3r  (6) 

3F  _ 3F  3X  3F  3Y 

30  3X  36  + 3Y  39 


The  transformation  functions  X and  Y for  any  given  problem  are  constructed 
analytically  to  provide  the  desired  non-uniform  mesh  spacing  in  the  physical 
plane.  That  is,  nodal  points  are  heavily  distributed  in  regions  of  large 
gradients  with  wider  spacing  elsewhere.  In  general,  X and  Y are  chosen  such 
that  the  magnitude  of  third  derivatives  of  the  dependent  variables  in  the 
computational  plane  does  not  become  large.  Since  X and  Y are  analytic  func- 
tions, their  derivatives  which  appear  as  coefficients  in  the  transformed 
equations  of  motion  via  Equations  (6)  are  also  analytic,  thereby  eliminating 
any  approximation  error  due  to  the  transformation  of  coordinates. 

At  each  interior  point,  (X,Y)  of  the  computational  plane,  the  dependent 
flow  field  variables  P,  R,  U and  V are  updated  at  each  time  step  by  means  of  a 
second-order  accurate  numerical  scheme  developed  by  MacCormack  (Reference  2). 
This  scheme  is  of  the  predictor-corrector  type  with  the  allowable  time  step- 
size  for  stable  computations  determined  by  the  Courant-Frieder ichs-Lewy  cri- 
terion (Reference  3). 

Construction  details  for  the  functions  X and  Y can  be  illustrated  quite 
simply  for  the  case  of  a symmetric  Joukowsky  airfoil  section  of  arbitrary  thick- 
ness h.  Conceptually,  the  transformation  proceeds  in  several  steps.  First, 
the  airfoil  section  denoted  by  b(0)  and  the  surrounding  outer  perturbation  wave 
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s(j,t)  sketched  in  Figure  2 are  mapped  from  the  physical  (r,6)  plane  to  an 
intermediate  (0,$)  plane  by  means  of  the  Joukowsky  transformation 


z = r,  + r + 


where 


= cr  + l£l 

and  r. . = 0 for  the  symmetric  airfoil  case.  This  mapping  produces  a smooth 
image  of  the  outer  perturbation  wave  denoted  by  S ( <t> , t ) in  the  c-plane  as 
shown  in  Figure  3 and  maps  the  airfoil  onto  a circle  of  radius  c provided 
a = c + e^.  Many  of  the  necessary  resolution  features  for  computing  complex 
airfoil  flow  fields  are  incorporated  in  this  transformation.  Note  that  the  two 
overlapped  rows  of  nodal  points  AB  and  CD  behind  the  airfoil  in  the  physical 
z-plane  facilitate  satisfying  boundary  conditions  at  the  X = 0 and  X * 1 
boundaries  of  the  computational  plane. 

Next  a normalizing  function  Z is  defined  in  the  C-plane  by 

Z = (9) 

S - c 

The  mapping  to  the  computational  plane  is  completed  by  constructing  functions 
of  the  form 

X = X(o,<f,a  ) 

Y = YCZ.Bj^)  °0) 

where  n^(r,9,t)  and  B^(r,0,t)  (i  * 1,2,...)  are  arbitrary  functions  chosen 
to  produce  a grid  network  in  the  physical  plane  having  the  desired  nodal  point 
distribution.  That  is,  points  should  be  heavily  concentrated  near  the  body 
and  its  wake  where  viscous  gradients  are  large,  near  the  outer  perturbation 
wave  where  pressure  gradients  may  be  steep,  and  in  the  vicinity  of  dividing 
streamlines  when  flow  separation  occurs  for  good  resolution  of  viscous  shear 
layers.  Moretti  (Reference  4)  has  successfully  computed  laminar  boundary  layer 
development  for  supersonic  blunt  body  flows  with  fewer  than  ten  points  dis- 
tributed across  the  layer  using  these  basic  concepts.  The  number  of  points 
required  to  accurately  predict  viscous  layer  development  will  be  studied  and 
determined  by  making  comparisons  of  the  TDC  program  results  with  the  results 
of  other  viscous  layer  calculations  such  as  Cebeci  (Reference  5). 
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as 


A simple  transformation  example  is  given  in  Reference  1 for  the  case  of 
inviscid  flow  past  a circular  cylinder.  The  normalizing  function  is  defined 


where  b is  the  cylinder  radius.  The  intermediate  mapping  to  the  -plane  is 
not  required  because  of  the  simple  geometry  (i.e.,  a = 0).  The  coordinate 
transformations  are  chosen  as 


and 


X = 


0 

TT 


Y 


1/2 


1 + (2Z  - 1)  tanh(q/2)  *1  ( 
1 - (2Z  - l)tanh(a/2) J j 


where  a(t)  is  arbitrary.  lines  of  constant  X are  simply  rays  in  the  physical 
plane.  Initially,  a is  chosen  such  that  an  approximately  linear  variation 
between  r and  Y along  each  ray  results.  As  the  physical  plane  expands  with 
time,  a linear  variation  is  not  sufficient  to  satisfactorily  resolve  the 
pressure  gradients  near  the  body  and  outer  perturbation  wave  without  using  an 
extremely  large  number  of  points.  By  allowing  a to  increase  with  time  in  such 
a manner  that  a constant  minimum  spacing  is  maintained  between  the  first  interior 
point  and  the  body  along  the  forward  stagnation  streamline,  fewer  points  are 
needed.  The  variation  in  mesh  point  distribution  at  various  times  is  shown 
qualitatively  in  the  sketch  below. 


0 y 1 

Although  the  grid  systems  resulting  from  the  coordinate  transformations 
described  above  present  many  computational  logic  problems,  the  resulting  smaller 
number  of  grid  points  required  to  adequately  resolve  the  flow  field  offers  a 
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considerable  saving  in  computer  storage  requirements  and  actual  computation 
time. 

2 . 3 Boundary  Conditions 

2.3.1  Initial  Conditions  - Initially  the  body  is  at  rest  followed  by  a 
smooth  acceleration  during  0 < t < t to  the  steady  velocity  VT  M • Hence,  at 
time  zero  values  of  the  dependent  variables,  P,  R,  U and  V are  zero  at  all 
points  of  the  flow  field.  The  motion  of  the  body  (r,0)  origin  relative  to  the 
initial  location  is  described  by 

x = -Vy”  M sinuit 

o ® 

Vy" 

x = (cosut  - 1) 

o u> 

for 

0 < t < t = 71  / 2u) 

s 

where  t is  arbitrarily  chosen, 
s 

At  the  start  of  motion,  a disturbance  wave  propagates  outward  from  the 

Initial  body  surface  location  at  the  speed  of  sound.  Up  to  a time  t^,  which 

is  small  compared  to  the  acceleration  time  tg,  the  physical  plane  grid  geometry 

is  held  stationary  with  respect  to  the  moving  body  origin  with  the  outer  boundary 

location  coincident  with  the  location  this  wave  front  will  have  at  time  t . 

o 

This  is  done  to  provide  a finite  size  physical  space  with  which  to  initiate  the 
computations.  Beyond  time  t , when  the  initial  perturbation  wave  just  reaches 
the  initial  grid  outer  boundary,  the  outer  grid  boundary  moves  along  with  the 
wave  front.  This  scheme  provides  for  better  definition  of  the  possibly  strong 
gradients  which  may  occur  along  portions  of  the  wave  and  allows  exact  far-field 
boundary  conditions  to  be  imposed. 

At  the  start  of  the  body  acceleration  (t  < t ),  the  velocities  and  pressure 

disturbances  are  small  and  the  wave  speed  can  be  assumed  constant  and  equal  to 

the  non-dimensional  free  stream  speed  of  sound, Vy*  Points  on  the  wave  front 

at  time  t will  lie  at  a distanceY/Y  t on  the  outward  normals  from  points  on 
o o r 

the  initial  body  surface  location  (Figure  4). 

Given  a point  on  the  body  surface  at  t = 0 located  at  (bQ,f)  relative  to 

the  origin  of  motion  (Figure  4),  let  the  corresponding  point  on  the  outward  normal 

at  the  wave  front  at  time  t be  at  (s  ) relative  to  the  origin  of  motion 

o o s,o 

and  be  located  at  (s,f7  ) relative  to  the  coordinate  origin  at  time  t . The  angle 

s o 
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between  the  outward  normal  and  the  radial  direction  at  (b  , S ) is  denoted  by 

o 


For  a given  value  of  5,  the  point  (s^.O^  n)  is  calculated  from  the  body  geometry 


by 


2.2  2 r- 

s =b  + yt  + 2 b V<  t cosr 
o o o o o 


(11) 


and 


VTt 

sin  (6  -5)  = — sinB 

9,0  s 

o 


(12) 


where 


, v2,-1/2 


b’ 


sin.-  = - — cosB 

D 

O 


and  primes  denote  differentiation  with  respect  to  Likewise,  for  a given 
value  of 
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2 2 

s + x + 2 s x cos0 
o o o o s,o 


sin (0  -0  ) = 

s ,o  s 


x sin" 
o llO 


A Newton-Raphson  iteration  on  S for  a given  value  of  0g  yields  a solution  of 

Equations  (11)  through  (14)  for  the  unknowns  <5,  0 , s and  s. 

M 6 s,o  o 

An  additional  level  of  iteration  may  be  required  depending  on  the  nature 

of  the  transformation  linking  the  physical  (r,9)  plane  to  the  computational  (X,Y) 

plane.  That  is,  a given  point  (X,l)  corresponding  to  a point  on  the  outer 

boundary  may  be  related  to  the  point  (s,0g)  by  relationships  such  as  those  given 

by  Equations  (7)  through  (10). 

3 s 

The  derivative  — required  in  the  course  of  the  flow  field  computation  for 

i3  0 

0 < t < t may  be  obtained  from  the  above  relationships  as 


s [sin(0g  0-6g)  “ F cos(0g-6)] 

F sin(6  -6)  - cos(0  -0  ) 
s s,o  s 


F = ^ 


b b'  +Vrt  (b'  cost!  - b [•? ' sinB) 
o o o o o 


s [N/yt  ts'  cosB  + s cos(0  -'5)1 
o o o s,o  J 


Also  required  is  the  derivative  — which  is  zero  for  0 < t < t . 

at  - o 


The  preceding  analysis  applies  to  that  part  of  the  outer  boundary  ahead 

of  the  outward  normal  from  the  initial  trailing  edge  location,  (x  + d),  or  for 

0 > 0_  (Figure  4).  For  9 < 0 the  boundary  is  a circular  wave  of  radius  \/i"t 
1 T o 

centered  at  the  trailing  edge  location.  The  outer  boundary  geometrv  and  deriv- 
ative for  a given  value  of  9 is  determined  as 

s 


(x  + d)  cosB  + \yt^  - (x  + d)^  sin^4 

n c ” n n c 


- s (x  + d)  sinO 
o s 

i — 2 7 r 

\ Y t -(x  + d)  sin  0 

’ o o s 


2.3.2  Wave  Front  Conditions  - The  outer  perturbation  wave  (Y  - 1)  sur- 
rounding the  body  is  compressive  upstream  and  expansive  downstream  of  the  body 


The  free  stream  velocity  components  (U  ,V  ) can  be  resolved  into  wave 
oriented  components  by 

V = (U  cost)  + V sin»))cosh  - (U  sin*)  - V cos<i)sinh  (lb) 

£ kxj  i *•  or 

and 

V = -(U  cost)  + V sint))sinh  - (U  sind  - V cos*i)cos,  (17) 

n oo  cw  au  < -jj 

Since  tangential  velocity  across  the  wave  front  is  conserved  whether  or  not 
the  wave  has  become  a discontinuity 

V = V cosh  - U sinh  (18) 

n 

where  U and  V are  the  radial  and  angular  components  of  velocity  immediately 
behind  the  wave  front.  Likewise,  the  normal  velocity  U,  immediately  behind 
the  wave  front  is 


(24) 


L 


and 


(25) 


If  the  wave  Is  expansive,  = u,,  so  that 


U = V. 


P = P 


and 


w = \ -rr 


(26) 


Differentiating  Equations  (23),  (24)  and  (25)  with  respect  to  T 


,+l 


(y+l)u^  J 


tT 


2u, 


3u, 


•T  2 t-1  IT 

"l  2 


(27) 


and 


iR 
■ T 


4'r(t  + l)u2  ,u! 


4> 


iu. 


2 > + (y-1)>iJ 


2 it 


(>+l)u 


IT 
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Differentiating  Equations  (20)  and  (21)  and  making  use  of  Equations  (27) 


1 _ Zl 

'T  (T 


<W 


and 


>T 


"2  - :“s.  « _ [tri  2v  1 "r 

T 'T  ’T  L'“  11 


Jf 

4 


so  that 


2ui  jwl 

2 _ ,-l  L >T  if  J 


>T  2 >-l  I IT 

it  - — L 

1 2 


(28) 
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r f , . * 


1R  = Ay 
IT  _ 2 

(>+l)u^u2 


'll.  2(>  + u")  f , 

t,  _ 1_  3W  >-l 

•T  , ' 2 IT  y+1 

(i+Duj  L 


2y  _§ 
2 >T 

(i+DujJ 


The  derivative — — appearing  in  these  expressions  may  be  evaluated  from 
Equation  (16)  as 


- sin 6 


where,  from  Equations  (15), 


/dU~  3 \ /dV„  30  \ 

_cosB  [dY-  + v.  W/  + sin\dT  ~v-^) 

+ <oU„  + VJ  (sint'>j  + sine  Lose  - u,  — j 

/du»  ae\  i 1 

\dY  + v~  W/+  (oV*  - u-}  It  (sln  > 


+ (ov  - U ) — (sin.  ) 

30  *33 


(sind)  = 


3 do 

cos  e — 
3T 


3t_  = 1 [ 3_  /3s\  _ 3_s"l 

3T  s [ 3T  \30  ) ° 3tJ 


3 / 3s\  32s  , 32s  30 


3T  ^ 30 ) 3t  30  2 3T 


The  unknown  shock  acceleration  — appearing  in  Equations  (28),  (29)  and 
(30)  can  be  determined  from  the  inviscid  equations  of  motion  by  forming  a 
characteristic  relation  valid  immediately  behind  the  wave  front  (Reference  6). 
Since  dY  0 along  the  wave 


1 5 


L 2l.-a2L 


Making  use  of  these  expressions,  the  equations  of  motion,  in  transformed 
computational  space  coordinates,  are  (see  Equations  (1)  through  (6)) 


DX 

21  + 

DY 

ip+  r 

Dt 

IX 

Dt 

Y 

U 

j.  I . 

oc 

S]- 

s 

S 

DO 

DX 

21  + 

DY 

DR  DY 

Dt 

Dt 

DY  )r 

U 

+ I - 

YX 

iX  = o 
>x 

s 

S 

)0 

DX 

21  + 

DY 

2 

DU  V 

Dt 

.x  + 

Dt 

DY  s " 

dU 

dV 

|+dX_U 
1 t r <X 


' Or  <Y  -r  >X  J 


= cosO  — — + sin9  — — 
dT  dT 


^V  + DX:)V  DYJV  yV 

>T  Dt  >X  Dt  DY  s ^ 


/l  . , jY  jP  \ 

10  DX  " ° lr  1 Y ) 


= cosO 


where 

ox  _.x  + y?x  + v ax 

Dt  Dt  Dr  s DO 


£Y  . Jl( a . oV  . is) 

Dt  >r  \ it  / 

By  means  of  Equations  (15),  (19)  and  (22) 

f - -L  ± (u,  - w) 

Dt  cos:-  >r 


Differentiating  Equation  (19)  with  respect  to  T 
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3Ur 

3T 


o / 'U  lV  \ 

cosb  ^ 7t  - 0 vr y 


cos^B  ~ (oU  + V) 


Combining  Equations  (34)  and  (35)  according  to  this  relation  gives 

A 


JUf.  + dv  + _g 3y  21 

3T  Dt  *Y  cosB  Jr  3Y 


< »6  i 


where 


js6  £( 


dU 


A = cosB  I (cos9  + osinO)  — — + (sinQ  - o cosb) 

dT 


dV, 

dT 


- cos  6 V 


_3o  _ DX  / 3U  9V\  Q^P(jX  £ jx\ I 

n 3T  Dt  \ 3X  3X  J ^ 3X  \3r  ~ s - /J 


VV 


Equation  (32)  may  be  written 


3P  DY  3P 


-v  30 f 

-Y—  ^ = B 


3T  Dt  )Y  cos6  3r  3Y 


(37) 


where 


B - - 5* 
Dt 


JP  _ rjX3U.U.l3X  3vl 
‘X  Y [ 3r  3X  s s 30  3xJ 


Forming  a characteristic  equation  from  Equations  (36)  and  (37), 
)U. 


Y 


<T  a 3T  + cost?  3r  (Uf  W a) 


r 3p  Y / 3U  „ v\ 

[iY-  a COS6  (Jy‘°-) 


= B 


where  a is  the  local  speed  of  sound.  Eliminating  the  term  by 

means  of  Equation  (33)  and  employing  Equations  (28),  (29)  and  (30),  this 
characteristic  equation  reduces  to 


iV  3 

~ = D E — f + F 
*T  [ 3T  ! 


1 38) 


where 


D = 


2u 


1 


2y(y+u2)  4^2(u2  - a) 

" 2 - S'1  (y+1 )au2  (y+l)au2u2 


-1 


2u 

■ 

" 

F i _ i. 

ri  _ 

2, 

E 2 v— 1 a 

)+l 

, , s 2 

v V 

('i+l)u1_ 

4 , (u2  - a ) 


( r+l)au^u2 


K 
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1 


F = [12~  a j_l  Tap  + lui  sr]  + i 

cosd  tr  I 3Y  a OY  I a 


A - B 


and  is  obtained  from  Equation  (31). 


3W 


Having  determined  the  shock  acceleration  — at  a given  node  point,  the 

aT 


shock  velocity  is  updated  by  a predictor-corrector  scheme  using  the  first- 
order  Taylor  approximation 

1W 


W(TKM)  = W(T)  + AT 


IT 


With 


the  updated  wave  velocity  W(T+AT),  updated  values  of  P,  R and  U 


immediately  behind  the  wave  may  be  determined  from  Equations  (20)  through 


jW 


(25).  For  the  expansion  portion  of  the  wave,  — is  obtained  by  simply  differentia- 


ting Equation  (26)  giving 
iV. 


*T 


while 


W(T+AT)  = Vf (T+AT)  - 


U (T+AT)  = V (T+AT) 


and 


P = R = 0 

At  each  nodal  point  along  the  outer  wave,  the  physical  plane  boundary  posi- 
tion can  then  be  updated  by  means  of  the  second-order  Taylor  expansion 


2 

s ( 9 , t + AT ) = s (0  , t ) + AT  ^ | AT2 

3T  2 aT2 


(39) 


where 


3s 

<T 


3s  3_s  39 
3t  + 30  3T 


(40) 


and 


2 2 2 

> s = s . 9 s 38 

,t2  2 3t 30  3T 

IT  it 


3_2_s  ( 3_8  \2  , 3s  3^8 
302V1T/  39  3T2 


(41) 


dX 


Since  -j-z,  r 0 for  a fixed  nodal  point, 
dT 


_i_0  = _ [dX  3X  3s  ] / 3x\_1 

IT  ” it  3r  3T J \ 30/ 
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-J 

-A 

v'ViT'T 


Finally,  from  Equations  (18)  and  (19), 

U = U;  cosF  - V sing 
V = U.  slnB  + V,  cose 

from  which  updated  values  of  U and  V may  be  determined  using  updated  values  of 
sine  and  cosg  computed  from  Equations  (15),  (39)  and  (42)  . 
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2.3.3  Body  Surface  Conditions  - Boundary  conditions  at  the  body  surface 
may  be  determined  in  a manner  quite  similar  to  that  for  the  outer  wave  con- 
ditions. That  is,  a characteristics  relationship  reflecting  the  fact  that 
the  normal  velocity  vanishes  at  the  surface  is  used  to  compute  surface  pressure. 
Hence,  a local  normal-tangent ial  coordinate  system  denoted  by  ( ,n)  in  Figure  6 


« 

L 


FIGURE  6.  RELATION  BETWEEN  NORMAL-TANGENTIAL  COORDINATES  (£,  r?)  AND 
RADIAL-ANGULAR  COORDINATES  (r,  n)  AT  BODY  SURFACE 


is  the  most  efficient  system  with  which  to  carry  out  the  computations.  The 
orientation  of  this  system  is  related  to  the  body  geometry  by 


cos  8 = (1  + a2)  1/2 


sin8  = - o(l  + o2)  1/2 


(43) 


where 


I 

b dO 


Also,  since  dY  --  0 along  the  body  surface 


1 3Y  = _ 0 3Y 
b 30  Dr 


(44) 


2. 3. 3.1  Inviscid  Flow  - Making  use  of  Equations  (43)  and  (44),  the 
normal  and  tangential  velocities  at  the  body  surface  measured  relative  to  the 
(r,0)  origin  fixed  in  the  body  are 
v = cosB  (U  - aV)  = 0 

and  (45) 

u = cosB  (oU  + V) 

Likewise,  the  inviscid  equations  of  motion,  in  transformed  computational  space 
coordinates,  are 


DP  A DX  DP  ^ 
dT  Dt  DX  + ' 


[ dy/du  iv\  3x  au  u . l ax  avl  = 0 

I Dr  ySY  ‘ ' DY  / 3 r DX  b b Dfi  3X  J 


(46) 
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_3U  U _1  _3X  _3V 

3X  b b 36  3X 


3R  DX  _3R  _3Y  /DU  3v\  JjX 

3T  Dt  3X  3r  ^3Y  ° 3Y ) dr 

jau  dx  au  v2  /ax  ap  dy  sp\  du-  , , dv~  . 0 

3T  Dt  3X  b + Q^ar  3X  + 3r  3Y  J ' dT  COSJ  + dT  Sln  ' 

d V HU 

i ax  bp  „ aY  ^p\  _ . 

b 39  3X  ~ 3r  3Y j dT  COSJ  " dT  Sln' 


(47) 


(48) 


(49) 


where 


DX 


3X  V 3X 

Dt  U 3r  b 36 


Differentiating  the  above  tangential  velocity  expression  with  respect  to  T 
and  using  the  fact  that  U = oV  on  the  body  surface 


3u 


3T 


= cosg 


/ 3U  3V  \ _ 1 

\°  3T  4 3T  / cos 


ay 

cosg  3T 


Combining  Equations  (48)  and  (49)  according  to  this  relation  gives  a tangential 
momentum  equation  of  the  form 


o dU  dV^ 

J V 2 00  GO 

— = cos  3 (o  cosO  - sinG)  + (o  sinG  + cosG) 

a I a I cl  I 


DX  / 3U  3V  \ _ / 3X1  3x\ 

Dt  \°  ax  4 ax  / Q\°  3r  4 b ae/ 


3P 

3X 


(50) 


w ith 


U = V 

Differentiating  the  normal  velocity  relation  with  respect  to  T gives 


av  al  au 

3T  = COs6l  3T  ' ° 


/ au  3v  \ _ 

^ DT  ° 3T  / ° 


Combining  Equations  (48)  and  (49)  according  to  this  relation  gives  a normal 
momentum  equation  of  the  form 


-y  + -4-  i*  » = cosB 

DT  cosg  3r  3Y 


dUa 

dV^ 

dT 

(cosO 

+ 0 

sin6)  ■+ 

dT 

( s i n 6 

“ 

O COS0) 

V2  , 

DX  | 

/ 3U 

- 0 

/ 3X 

0 

ax  \_ap  1 

cos^B 

Dt  1 

\ DX 

3X  ) 

X 

V 3r  - 

b 

36  / ax 
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Forming  a characteristic  equation  with  this  result  and  Equation  (Ah), 


jY 

' IE  _ 

y(2U  _ 

, iv  V 

DX 

3P 

or 

cos 6 3Y 

r\3Y 

*Y  / 

Dt 

3X  ” 1 

+ U + 

1) 


i ax  av 

b 30  ax 


(cos0  + o sin0) 


dV-  . . V2  DX  / 3U 

+ d'T  (sin  ' ' cos  ) + bVo^B  - nt  ( VX 


/ 2X  o 1X\  3P 
^3r  b 39 ) ax 


where  a is  the  local  sonic  velocity.  Finally,  Equations  (A6)  and  (47)  may 
be  combined  to  give 


1 

'lE  + 

DX 

(lE 

3R  \ 

Y 

[_3T 

Dt 

^ 3X 

r 3X  J 

which  reflects  the  fact  that  entropy  is  convected  along  streamlines. 

Updated  values  of  P,  R,  U and  V at  points  on  the  body  surface  can  be  deter- 
mined by  means  of  Equations  (50),  (51)  and  (52).  However,  special  consideration 
must  be  given  to  points  of  surface  discontinuity  such  as  sharp  trailing  edges. 
This  will  be  illustrated  below  for  the  case  of  a non-lifting,  symmetric  Jnukowsky 
airfoil. 

The  Joukowsky  transformation  defined  by  Equations  (7)  and  (8)  gives  rise 
to  a geometric  singularity  at  the  cusped  trailing  edge  which  in  turn  induces 
singular  behavior  in  the  spatial  derivatives  of  the  flow  variables  P,  R,  U and  V. 
This  is  evident  by  expanding  in  terms  of  body  surface  polar  angle  0 near  the 
trailing  edge  the  quantity  a defined  by  Equations  (43).  That  is, 


where 


o - Kl  <r1/3  + k2  e1/3 


k2  --!  KK.  (tt t) 


+ K3  0 + • 


k = ( c + e 


."W 


Since,  from  Equations  (50),  U = aV  along  the  body  surface,  in  order  that  the 
trailing  edge  velocity  remain  finite,  V must  possess  an  expansion  of  the  form 
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V = Ax  0 ' + A,  0 + A3  I + 


(5 


so  that 


U = AlKj  + (A3K7  + A.jICj  ) 2/t  + CA?  K J + A,K,  + A^)  + •••  (5 


I'h  is  implies 


U- 


'TK 
A , ~ - 

1 Ki 


where  UyE  denotes  the  trailing  edge  velocity.  Analogous  expansions  may  be 
assumed  for  P and  R having  the  form 


P + PTp  + Bj  02/3  + B.,  04/J  + 


(5 


and 
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Along  the  body  surface 
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so  that  the  inviscid  equations  of  motion  may  be  written 


V 

dP 

4. 

' m U 

1 3V 

b 

dO 

> 

y\ 

3r  b 

b 36 

V 

dR 

4. 

iy 

+ u + 1 

3 V 

b 

de 

r 

3r 

+ b + b 

30  ' 

V 

dU 

v2 

+ o 3P  _ 

dU 

CO 

b 

de 

b 

+ Q 3r  ~ 

dt  ' 

V 

dV 

1 

UV 

+ 2 3P  _ 

dU 

b 

de 

b 

b 30 

dt 

(; 


(f 


sin  > 


(' 


3V 


Since  — = 0 at  the  trailing  edge  for  the  non-lifting,  symmetrical  ui 
foil  section,  Equations  (53)  and  (bl)  imply 
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and  therefore,  by  Equations  (55)  and  (57), 
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Equations  (58),  (59)  and  (60)  are  not  amenable  to  finite  difference 
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computation  at  the  trailing  edge  since  the  derivatives  — , — - and  — are 
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singular  at  this  point  although  the  terms  V— , V—  and  Vjp  remain  finite. 


This  difficulty  can  be  circumvented  by  considering  instead  the  derivatives 
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obtained  from  Equations  (53)  through  (56).  The  problem  terms  appearing  in  the 
equations  of  motion  may  be  written 
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so  that  at  the  trailing  edge  where  6=0, 
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and  from  Equation  (62) 
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The  term  — - in  Equations  (58)  and  (59)  mav  be  written  in  computational 
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space  coordinates  as 
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where  — = 0 has  been  assumed  since  body  surface  points  are  being  considered. 
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Since  the  Joukowsky  transformation  defined  by  Equation  (7)  is  analytic,  the 
Cauchy-Riemann  conditions 
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must  be  satisfied.  Equation  (44)  is  equivalent  to 


since  p E c along  the  body  surface.  Therefore, 
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and,  from  Equations  (57),  (65),  (66)  and  (67), 
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Equation  (64)  may  now  be  written 
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Analysis  of  the  Joukowsky  transformation  yields  the  result,  valid  on  the 
airfoil  surface , 
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The  limiting  value  of  the  term  o — in  Equation  (69)  can  be  determined  by 
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L'Hospital's  rule.  That  is,  using  Equation  (68), 
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Therefore,  the  limiting  value  of  Equation  (69)  becomes 
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where  7--  = 0 has  been  assumed  because  of  symmetry. 
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In  a similar  manner,  the  term  - — in  Equations  (58)  and  (59)  may  be 
written 
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whe  re 
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By  means  of  Equations  (65)  and  (66), 
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Again  using  L'Hospital's  rule  (refer  to  Equation  (71)), 
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Therefore,  the  limiting  value  of  Equation  (74)  becomes 
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Making  use  of  Equations  (63),  (72)  and  (75),  the  equations  governing  the 
flow  properties  at  the  trailing  edge  are 
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All  terms  in  these  expressions  can  be  evaluated  at  the  trailing  edge  in  a 
straightforward  manner. 
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Since  MacCormack's  scheme  makes  use  of  both  forward  and  backward  spatial 
differencing,  it  is  consistent  to  develop  relations  analogous  to  the  above 
Equations  (76)  based  on  an  analysis  of  flow  conditions  immediately  downstream 
of  the  trailing  edge.  For  flow  along  the  trailing  edge  streamline,  the  equa- 
tions of  motion  reduce  to 
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The  singular  nature  of  the  trailing  edge  geometry  again  induces  singular 
behavior  in  the  derivatives  of  the  flow  variables  P,  R and  U appearing  in  these 
equations . 

From  the  J oiitowsky  transformation,  it  can  be  shown  that 
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along  the  trailing  edge  streamline,  where  b denotes  the  value  of  r at  the 
trailing  edge.  This  latter  result  suggests  expansions  of  the  form 
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Equations  (77)  and  (78)  may  be  combined  to  give 


<P  :,R  X 11 

— - v— - + U 
it  ' ,)t 


/ dP  JR\ 

fa  ■ Y 


If  this  result  is  to  be  well-behaved  at  the  trailing  edge,  then  the  above 
expansions  lead  to 
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likewise,  if  Equations  (77)  and  (79)  are  to  be  non-singular  at  the  trailing 


edge, 
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With  these  results,  the  governing  equations  at  the  trailing  edge  become 
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From  Equations  (81)  and  (83)  evaluated  at  r = b, 
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Considering  symmetry  and  Equations  (65)  and  (80), 


1 jV  = 1 ,V  _tX  = 1 Al  jX  _ ±V  dX  1 _3£ 
r r 'X  r 'X  • ; )0  'X  '•••  p t 


'X  iV  1 

£ j 

fr2  - h2  V1'2  2'r 

“ 

^ X r - 2- 

r 

b 1 

l b2  ) + - 

Relating  this  result  to  the  expansion  for  “ given  in  Equations  (81), 
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Using  Equations  (85),  (86)  and  (87),  the  equations  of  motion  at  the 
trailing  edge,  valid  for  downstream  differencing,  are 


— = i F - I U_  (P-P  ) + Y(U  - U )| 
at  lr  L TE  TE  TE  J 

Ti  - F - Y7  [l,TE  <“  - V + u - “te]  <' 

do,  '2 b qte2  /ixwvr,  (d 

dt  <.c2«4.  - ,qte>2  '*)  [ > v 

7f»TE  <"  - »Tr>  * «r,  C - P,*>1 


and 


>U 

at 


,3 


2. 3. 3. 2 Laminar  Viscous  Flow  - The  analysis  of  viscous  boundary  conditions 
is  simplified  somewhat  by  the  no-slip  condition.  If,  in  addition,  a constant 
temperature  surface  is  assumed  equal  to  free  stream  temperature,  the  boundary 
conditions  to  be  applied  to  Equations  (1)  through  (5)  become 

U = V = 0 

Q = 1 = eP_R  (89) 

or 

P = R 

This  last  condition  reduces  the  solid  boundary  analysis  to  that  of  predicting 
only  the  surface  pressure  which  can  be  accomplished  most  accurately  by  con- 
structing a characteristic  equation  from  Equations  (2),  (3)  and  (4). 

Following  a procedure  similar  to  that  outlined  for  inviscid  flow,  a 
normal  momentum  equation  is  formed  by  combining  Equations  (2)  and  (3)  accord- 
ing to  the  first  of  Equations  (45).  This  result  is  then  combined  with  Equa- 
tion (4)  to  give  the  characteristic  relationship 
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In  these  expressions,  use  has  been  made  of  Equation  (44)  and  the  additional 
surface  conditions 
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3X  ' 3X 


2 2 
3 r 3 R 


obtained  from  Equations  (89).  The  quantities  F,  G and  H appearing  in  Equations  (91) 
represent  the  viscous  terms  in  Equations  (2),  (3)  and  (4)  which  have  been 
assumed  to  act  as  forcing  terms  when  forming  the  characteristic  relation  of 
Equation  (90).  This  technique  has  also  been  used  in  Reference  7. 

Equation  (1)  may  be  written  for  constant  surface  temperature  as 
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which  can  be  incorporated  into  Equation  (90)  to  eliminate  the  term 
giving 
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This  result  may  then  be  used  to  update  surface  pressure. 

The  special  treatment  required  at  points  of  surface  discontinuity  will 
again  be  illustrated  for  the  cusped  trailing  edge  of  a non-lifting  Joukowsky 
airfoil.  For  surface  points  near  the  trailing  edge.  Equations  (1)  and  (4) 
may  be  combined  in  computational  space  variables  as 
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By  means  of  Equations  (65)  and  (66)  and  the  fact  that  r-7  and  — - both 
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Differentiation  of  the  Cauchy- Riemann  conditions  gives 
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so  that  Equation  (94)  becomes 


In  order  that  this  term  be  non-singular  at  the  trailing  edge,  the  quantity 
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must  vanish.  Equation  (93)  evaluated  at  the  trailing  edge  then  becomes 
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Double  application  of  L' Hospital's  rule  and  Equation  (68)  gives 
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Therefore,  the  governing  equation  at  the  trailing  edge  based  on  upstream 
influence  is 
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Along  the  trailing  edge  streamline  the  continuity  and  energy  equations 
simplify  to 
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Considering  the  term — in  the  above  equations  and  employing  Equations 
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(65)  and  the  fact  that  the  flow  and  coordinate  geometry  are  symmetric 
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From  the  inviscid  analysis,  the  derivative  ~ given  by  Equations  (80)  is  singular 
at  the  trailing  edge.  Applying  I.’ Hospi tal 's  rule  to  the  above  expression  gives 
the  trailing  edge  value 
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which  is  well-behaved  and  differs  from  the  inviscid  result.  This  implies  that 
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This  result  may  be  used  to  simplify  Equation  (96).  By  similar  arguments,  the 
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Considering  the  remaining  terms  in  Equations  (98)  and  using 
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must  vanish  at  the  trailing  edge.  This  requirement  is  ensured  by  Equation  (100) 
and  its  preceding  arguments. 

Using  Equation  (99),  Equation  (101)  evaluated  at  the  trailing  edge  then 
becomes 
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Application  of  L'Hospital's  rule  gives 
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Therefore,  from  Equations  (97)  the  governing  equation  at  the  trailing  edge 
based  on  downstream  influence  is 
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2.3.4  Symmetry  or  Identity  Conditions  - Boundary  conditions  at  X = 0 
and  X = 1 are  much  simpler  to  determine.  If  the  body  and  flow  field  are 
symmetrical,  an  extra  row  of  nodal  points  is  added  below  the  symmetry  boundary 
both  fore  and  aft.  Symmetry  or  anti-symmetry  conditions  on  the  flow  variables 
are  then  imposed  along  both  of  these  extra  rows  of  points.  If  the  flow  field 
is  not  symmetric,  the  X = 0 and  X = 1 rows  of  nodal  points  are  made  to  overlap 
by  one  grid  interval  at  some  point  in  the  field  (see  Figure  2).  Identity  con- 
ditions are  then  imposed  on  the  two  overlapping  rows  of  points. 
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2 . 4 TDC  Computer  Program  Organization 

A flow  chart  of  the  TDC  program  is  presented  in  Figure  7.  Subroutine 
INITAL  handles  the  input  and  storage  of  the  body  geometry  and  its  surface  de- 
rivatives and  sets  the  initial  values  (normally  zero)  of  the  dependent  vari- 
ables P,  R,  U and  V.  For  some  small  time  tQ,  the  outer  wave  location  is 
computed  (see  Section  2.3.1)  by  assuming  a propagation  velocity  equal  to  free 
stream  sonic  velocity.  This  wave  location  initially  serves  as  the  outer 
boundary  of  the  physical  plane  and  moves  with  the  body  for  0 < t <_  tQ.  For 
t > tQ,  the  physical  plane  is  allowed  to  expand  in  a normal  fashion  with  the 
outer  wave  location  and  properties  determined  by  subroutine  OUTERW  as  explained 
earlier  (see  Section  2.3.2.). 

Subroutine  DIFFER  contains  the  finite  difference  equivalent  of  Equations 
1 through  5 and  updates  the  dependent  variables,  P,  R,  U and  V at  each 
interior  point  of  the  computational  space.  The  allowable  time  step-size  is 
computed  in  this  subroutine  by  means  of  the  Courant-Friederichs-Lewy  (CFL) 
rule  as  modified  for  viscous  flow  by  Moretti  (Reference  4).  That  is, 


AT  = 


kdmin 


yJv2  + V2  + a + 

min 


(103) 


where  dm^n  is  the  miiumum  spatial  increment  in  the  flow  field,  a is  the  local 
speed  of  sound,  and  k is  a safety  factor  less  than  unity. 

Subroutine  TRANSF,  called  by  DIFFER,  analytically  computes  the  derivatives 
of  the  transformation  functions  X and  Y employed  in  mapping  the  physical  space 
onto  the  computational  space.  The  choice  of  these  transformation  functions  is 
entirely  arbitrary,  being  constructed  to  suit  the  point  distribution  require- 
ments of  each  given  type  of  problem.  Subroutine  BNDARY  computes  conditions  at 
the  boundaries  of  the  computational  space  shown  in  Figure  1.  The  nature  of 
these  boundary  conditions  will  vary  somewhat  with  each  given  class  of  problems 
and  the  assumptions  concerning  the  type  of  flow,  such  as  inviscid  flow,  iso- 
thermal or  adiabatic  wall,  etc. 

Subroutine  OUTPUT  writes  out  the  dependent  variables,  outer  wave  geomet- 
rical characteristics,  and  other  parameters  of  interest  after  each  chosen  in- 
crement in  time  steps.  All  of  the  information  in  machine  storage  may  also  be 
stored  on  tape  after  an  arbitrary  number  of  time  steps  so  that  the  computations 
may  be  resumed  at  a later  date. 
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FIGURE  7. 


ORGANIZATION  OF  TIME  DEPENDENT  COMPUTATION  (TDC) 
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3.  EXTENSION  OF  TDC  PROGRAM  FORMULATION  TO  TURBULENT  FLOW 


The  Time  Dependent  Computation  (TDC)  computer  program  formulation  for 
inviscid  and  viscous  laminar  flow  is  given  in  Section  2.  The  turbulence  model 
selected  for  incorporation  into  the  program  is  fashioned  from  the  Extended 
Energy-Dissipation  model  of  Spalding  (Reference  8)  with  additional  corrections 
for  compressibility  as  developed  by  Wilcox  for  the  Saffman  model  (Reference  9). 
3.1  Spalding  k-e  Turbulence  Model 

The  turbulence  model  selected  employs  an  equation  for  both  the  turbulence 
energy  k and  the  turbulence  energy  decay  rate  e.  The  introduction  of  the  decay 
rate  eliminates  the  necessity  of  prescribing  a length  scale,  so  that  the  turbu- 
lent eddy  viscosity  is  given  by 


C p - 
U e 


(104) 


where  p is  the  density  and  C^  is  the  turbulent  viscosity  coefficient.  From 
Reference  9,  the  equations  in  rectangular  coordinates  governing  k and  c are 


(105) 


(106) 


where  u is  the  laminar  viscosity,  <|>  is  the  laminar  dissipation,  u^  represents 

velocity  components  and  o,  , o , C , C . and  f;  are  empirical  constants.  Com- 

Keeicz  3u£ 

pressibility  effects  are  accounted  for  by  the  term  £pk  - — in  Equation  (105). 

3 X . 


The  use  of  in  the  production  terms  of  Equations  (1^5)  and  (106)  is  a 
generalization  from  the  boundary  layer  approximation  to  two-dimensional  flow 
of  the  corresponding  terms  in  Spalding's  model  and  are  identical  to  the  terms 
in  Wilcox's  compressible  form  of  Saffman' s model  (Reference  10).  The  turbulent 
viscosity  coefficient  in  Equation  (104)  is  evaluated  from 


^ ' -OT  f(i) 

where  f is  given  in  Figure  1 of  Reference  8.  The  argument  of  this  function  Is 


7-/\  <*]  1 

where  t is  the  turbulent  shear  stress  and  the  integrations  are  carried  out 
across  the  viscous  region. 
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The  values  of  the  empirical  constants  appearing  in  Equations  (105)  and  (106) 
are  given  in  Reference  8 as 


o.  = 1.0 
k 

Ccl  = 

1.40 

o =1.3 

e 

Cc2- 

1.94 

The  value^of  £ was  determined  in  Reference  10  as 
5 « 2.5 

Hoffman  (Reference  11)  has  examined  the  behavior  of  the  Spalding  k-e  turbulence 
model  for  consistency  at  a turbulent/non-turbulent  interface  and  recommends  the 


: different 

set  of  constants 

o,  = 2.0 

C . = 1.81 

k 

cl 

o = 3.0 

e 

C , = 2.00 
e2 

Improved  agreement  with  experimental  results  for  fully  developed  channel  flow  was 
obtained  with  these  values. 

3. 2 Equations  of  Motion 

The  turbulence  model  described  above  introduces  additional  terms  and  modi- 
fications into  the  Navier-Stokes  equations  as  presented  in  Section  2 to  make 
them  valid  for  turbulent  flow.  The  dependent  variables  P,  R,  U and  V now  denote 
mean  instead  of  instantaneous  values  of  pressure,  density  and  velocity. 

Under  the  assumptions  of  a perfect  gas  with  constant  laminar  viscosity  and 
thermal  conductivity,  the  equations  of  motion  in  polar  coordinates  are 


dt  9r  r 96  1 

an  au  v au  c^au  „ 3P 

3t  3r  r 39  r dt  dt  3r 

♦ )♦(««£ 


(107) 


(108) 
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and 
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The  symbols  used  in  these  equations  are  defined  as  follows: 

E - Turbulence  energy  decay  rate  non-dimensionalized  by 
K - Turbulence  energy  non-dimensionalized  by  P /d 

OO  CO 

L - Reference  length 


*(£) 


3/2 


P-R 


P - Log  (p/pj 

Q - Non-dimensional  temperature,  e‘ 

R - Log  (p/o^) 

U - Radial  velocity  non-dimensionalized  by  Vp/p 


V - Angular  velocity  non-dimensionalized  by  \/p ,„/ P,, 
p - Static  pressure  (mean) 

r - Radial  coordinate  non-dimensionalized  bv  L 


(mean) 

(mean) 


t - Time  non-dimensionalized  by  L Vr>  / i> 
<t>  - Laminar  dissipation  function 
p - Mass  density  (mean) 

0 - Polar  angle 
Y - Ratio  of  specific  heats 
5 - Turbulence  model  constant 
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k 


M - Free  stream  Mach  number 

ao 

Re  - Free  stream  Reynolds  number 
Pr^  - Laminar  Prandtl  number 

Pr^  - Turbulent  Prandtl  number 

- Longitudinal  free  stream  velocity  component 

- Vertical  free  stream  velocity  component 

P^  - Free  stream  static  pressure 
p - Free  stream  mass  density 

- Turbulent  eddy  viscosity  non-dimensionalized  by  L 

- Turbulent  viscosity  coefficient 

a.  - Turbulence  model  constant 

k 

a - Turbulence  model  constant 

e 

C , - Turbulence  model  constant 

El 

C „ - Turbulence  model  constant 

e2 

The  finite  difference  equivalent  of  the  above  equations  of  motion  are 
contained  in  subroutine  DIFFER  (see  Section  2.4)  where  the  dependent  variables 
are  updated  at  interior  points  of  the  computational  space.  If  conditions  are 
such  that  at  some  point  in  the  flow  field  transition  to  turbulent  flow  occurs, 
the  above  equations  are  employed  to  determine  flow  properties  at  all  downstream 
and  all  transverse  locations  likely  to  be  affected. 

3.3  Boundary  Conditions 

Turbulent  flow  boundary  point  calculations,  to  a certain  extent,  are  the 
same  as  for  laminar  flow  provided  the  turbulent  flow  is  confined  to  the  body 
boundary  layer  and  wake  region.  Initially,  the  dependent  variables  P,  R,  U, 

V,  E and  K are  each  zero  with  the  outer  wave  front  boundary  location  determined 
as  described  in  Section  2.3.1.  Properties  immediately  behind  the  moving  outer 
perturbation  wave  can  again  be  determined  by  assuming  inviscid  conditions  as  in 
Section  2.3.2. 

It  is  advantageous  to  avoid  calculating  turbulence  properties  directly  very 
near  a solid  boundary  using  Equations  (111)  and  (112)  because  of  the  extremely 
small  grid  spacings  required  to  obtain  adequate  resolution  of  the  rapid  varia- 
tion of  k and  e across  the  laminar  sublayer.  This  difficulty  can  be  avoided  by 
imposing  fictitious  values  of  k,  e and  slip  velocity  at  the  surface  such  that 
reasonably  accurate  spatial  derivatives  normal  to  the  surface  can  still  be 
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obtained  by  conventional  finite  differencing  with  not  too  fine  a mesh.  The 
necessary  approximations  to  the  turbulent  flow  properties  near  the  surface 
can  be  obtained  using  "law-of-the-wall"  relationships.  These  relations  are 


r 


+ 1 i . +x 

U = A 108  ^By  ^ 


and 

+ 1 

e + 

Ay 

where  by  definition 


u+  b *- 
u 

T 


A = .4187 
B = 9.793 


(113) 

(114) 

(115) 

(116) 

(117) 

(118) 

(119) 

(120) 


and  y is  the  distance  normal  to  the  surface,  u is  velocity  parallel  to  the 


surface  and  p 


and  denote  surface  values  of  density,  laminar  viscosity 


and  shear  stress. 

Because  of  the  logarithmic  variation  of  u and  the  inverse  power  variation 
of  c,  the  above  equations  are  not  valid  at  the  surface.  However,  several 
simplifying  assumptions  can  be  made  to  achieve  the  desired  approximate  result. 
First,  the  shear  stress  t is  assumed  to  be  constant  near  the  surface.  Next, 
values  of  u+,  k+  and  e+  are  determined  at  a point  p (see  Figure  8)  near  the 
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inner  edge  of  the  " law-of- the-wa 1 1"  region  and  surface  values  are  extrapolated 
linearly  from  this  point.  The  point  p is  chosen  to  coincide  with  the  location 
at  which  the  velocity  variation  of  Equation  (113)  matches  the  laminar  sublayer 
law  u+~  y+ , namely  y+  = 11.63.  The  first  grid  point  outward  from  the  surface 

P + 

should  lie  near  the  outer  limit  of  the  " law-of- the-wall"  region  where  y is 
approximately  100  to  200  and  the  wall  shear  stress  will  be  assumed  to  be  the 
same  as  at  this  point. 

The  necessary  slopes  for  extrapolating  to  the  fictitious  surface  values 
can  be  obtained  by  differentiating  Equations  (113)  through  (115).  The 
surface  values  of  u+,  k+  and  e+  are  then  given  by 


u+  = I log  (By+)  -1  = 8.921 

w A p 


k+  = c-1/2 

w Vi 


3.333  f 


+ 2 

e = — = .4107 

w •+■ 

% 


From  Equations  (116)  through  (120),  the  surface  values  of  k,  c and  velocity 
in  non-dimensional  form  are 


K = 3.333  f ' t 

w w 
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f t e w 
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U = aV 
w w 


1 db 

11  b dO 


where  R denotes  the  surface  value  of  R and  b is  the  radial  coordinate 
w 

describing  the  body  surface.  Note  that  the  slip  velocity  is  non-zero  for 
the  turbulent  boundary  layer  only. 

Surface  values  of  pressure  for  turbulent  flow  can  be  approximated  using 

the  characteristics  relationship  developed  in  Section  2. 3. 3. 2.  The  assumption 

underlying  such  a prediction  is  that  the  turbulence  quantities  in  Equations 

(108),  (109)  and  (110)  are  negligible  at  the  surface.  The  analysis  of 

Section  2. 3. 3. 2 assumes  a constant  wall  temperature  equal  to  free  stream 

temperature  which  implies  equality  of  P and  R along  the  surface.  At  a sharp 

trailing  edge,  the  above  approximation  using  the  "law-of-the-wall"  relations 

will  not  be  valid.  However,  values  of  E , K , U and  V can  be  approximated 

WWW  w 

by  simple  extrapolation. 
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4.  INVISCID  SAMPLE  CASES 

Several  exact  inviscid  two-dimensional  flow  field  methods  are  available 


for  the  purpose  of  validating  the  results  of  the  TDC  Program: 

° Douglas  Neuman  Program  (Reference  12)  for  incompressible 
potential  flow 

° Garabedian-Korn  Analysis  Program  (Reference  13)  for 
subcritical  or  supercritical  airfoil  flow 
° Catherall-Sells  Program  (Reference  14)  for  subcritical  flow. 

In  this  section,  calculations  by  the  Catherall-Sells  Program  are  used  as 
the  basis  for  comparison.  This  method  uses  a relaxation  technique  to  solve 
the  subsonic,  inviscid  flow  equations  in  a computational  space  consisting  of 
a polar  grid  inside  a circle.  Comparisons  of  TDC  Program  calculations  with 
the  steady-state  calculations  obtained  from  the  Catherall-Sells  Program  were 
used  to  determine  the  number  of  grid  points  needed  for  accurate  results,  the 
time  of  convergence  to  steady  state  values  at  the  body  surface,  and  the  rate 
of  growth  of  the  steady  state  region  around  the  body. 

4. 1 Non-Lifting  Flow  PasL  a Circular  Cylinder 

Sample  results  for  inviscid  flow  past  a circular  cylinder  of  unit  radius 
are  presented  in  Figures  10  and  15  for  several  subcritical  free  stream  Mach 
numbers.  The  grid  system  for  these  cases  consisted  of  rays  of  points  from  the 
center  of  the  cylinder  equally  spaced  over  the  upper  surface  in  the  angular 
direction  as  shown  in  Figure  9.  The  radial  distribution  of  points  along  each 
ray  is  such  that  points  are  concentrated  near  the  body  and  near  the  outer 
perturbation  wave  where  gradients  tend  to  be  largest.  The  polar  angle  6 is 
measured  from  the  downstream  axis  while  r is  measured  outward  from  the  center 
of  the  cylinder.  The  pressure  coefficient  is  defined  as 
c = iQR  P - 1 
P 1/2  Y M2 

QD 

The  transformation  functions  X and  Y for  this  case  were  chosen  similar 
to  those  of  Reference  1,  namely 


X = - 
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Y = 2 + 2i'°S 
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FIGURE  9.  GRID  GEOMETRY  FOR  CYLINDER  FLOW 


Catherall  Sells  Solution 
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FIGURE  10.  CYLINDER:  SURFACE  PRESSURE  DISTRIBUTION 


Catherall  Sells  Solution 


Incompressible  Solution 


FIGURE  11.  CYLINDER:  NEAR  FIELD  PRESSURE  DISTRIBUTION 


The  derivatives  of  these  functions  required  in  the  course  of  the  numerical 
computation  are 
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The  parameter  a is  initially  chosen  to  give  an  approximately  linear 
variation  between  r and  Y (e.g.  , a = 1.0).  As  the  physical  space  expands 
with  time,  a is  adjusted  according  to 

dmin  1 ) + tanh[n  (AY  _l_/_2)] 

s(ir,t)  - 1 2 | tanh(a/2) 
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La  maintain  a constant  spacing  d between  the  first  interior  point  and  the 

'tun  r 

body  along  the  forward  stagnation  streamline.  The  time  derivative  of  j re- 
quired in  the  above  expressions  is  obtained  from  this  equation  as 


da 

dt 
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Figure  10  compares  the  calculated  surface  pressure  distribution  for  free 
stream  Mach  number  of  0.1  with  the  exact  incompressible  solution  and  the 
Catherali-Sells  solution.  Reasonable  accuracy  is  obtained  with  only  21  upper 
surface  points  and  31  points  in  the  radial  direction.  The  non-dimensional 
time  T = 12.8  ccrrsnonds  to  1000  time  stens.  Figure  11  shows  the  calcu- 
lated pressure  distribution  near  the  cylinder  compared  with  the  Catherali- 
Sells  solution  and  indicates  a steady  flow  region  extending  outward  approxi- 
mately seven  radii.  The  effect  of  the  coordinate  stretching  function  Y is 
clearly  e”ident  in  that  the  point  spacing  near  the  body  where  gradients  art- 
large  is  much  smaller  than  farther  out  in  the  flow  field. 

Figure  12  compares  the  calculated  surface  pressure  distribution  for  free 
stream  Mach  number  of  0.2  with  the  Catherali-Sells  solution.  The  incompres- 
sible solution  is  also  shown  for  reference.  Excellent  agreement  (.less  than 
0.5%  error) i^  achieved  in  this  case  with  41  points  in  the  angular  direction 
and  51  points  in  the  radial  direction.  Figure  13  shows  the  near-field  pres- 
sure distribution  compared  with  the  Catherali-Sells  solution  and  indicates  a 
steady  flow  region  extending  outward  approximately  eight  radii.  Figure  14 
presents  a time  history  of  the  upstream  and  downstream  stagnation  pressures 
and  indicates  at  which  time  steady  conditions  first  begin  to  appear  near  the 
body.  Figure  15  compares  the  calculated  surface  pressure  distribution  for 
free  stream  Mach  number  of  0.3  with  the  Catherali-Sells  solution,  again  with 
good  agreement.  Critical  free  stream  Mach  number  for  a circular  cylinder  is 
approximately  0.4. 

4. 2 Lifting  Flow  Past  a Circular  Cylinder 

Sample  results  for  inviscid  flow  past  a circular  cylinder  with  arbitrary 
circulation  are  presenteu  in  Figures  1(>  through  21  for  suberiticai  free  stream. 
Mach  numbers  of  0.1  and  0.2.  The  grid  arrangement  and  transformation  func- 
tions X and  Y were  similar  to  those  for  the  above  non-lifting  case  except  for 
extending  around  the  lower  half  of  the  cylinder  also,  as  shown  in  Figure  9. 
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FIGURE  17.  CIRCULAR  CYLINDER  WITH  CIRCULATION  (INVISCID) 
NEAR-FIELD  PRESSURE  DISTRIBUTION 
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FIGURE  20.  CIRCULAR  CYLINDER  WITH  CIRCULATION  (INVISCID) 
NEAR-FIELD  PRESSURE  DISTRIBUTION 
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Due  to  the  absence  of  any  flow  symmetry,  the  grid  system  has  an  overlap  region 
near  0=0  for  simple  specification  of  boundary  conditions  at  X = 0 and  X = 1 
of  the  computational  space. 

The  TDC  solutions  were  obtained  by  assuming  that  initially  the  cylinder 
was  surrounded  by  a compressible  vortex  flow  of  arbitrary  circulation  concentric 
with  the  cylinder.  The  cylinder  was  subsequently  accelerated  smoothly  from 
rest  to  free  stream  velocity  with  the  vortex  center  coincident  with  the  center 
of  the  cylinder.  Outside  of  the  initial  perturbation  wave,  the  combination 
of  vortex  and  rectilinear  flow  remained  unaffected  by  the  motion  of  the  cylinder 
so  that  correct  far-field  boundary  conditions  could  be  simply  imposed.  Care- 
ful, systematic  numerical  studies  have  indicated  that  this  is  the  physically 
realistic  manner  in  which  to  impose  lift  in  the  development  of  a time  dependent 
inviscid  flow. 

The  surface  pressure  distributions  shown  in  Figures  16  and  19  were  ob- 
tained with  40  points  distributed  uniformly  around  the  cylinder  and  31  points 
in  the  radial  direction.  Although  the  point  distribution  density  was  approxi- 
mately half  that  employed  in  Figures  12  and  13 , good  agreement  (approximately 
1%  error)  with  the  Catherall-Sells  solution  was  still  obtained.  The  value  of 
circulation  around  the  cylinder  was  chosen  to  fix  the  steady  state  rear  stag- 
nation point  at  the  preselected  angular  position  designated  by  6^.  Figures  17 
and  20  show  the  near-field  pressure  distributions  compared  with  the  Catherall- 
Sells  solution  and  indicate  a steady  flow  region  extending  outward  approxi- 
mately six  radii.  Time  histories  of  lift  and  drag  coefficient  are  shown  in 
Figures  18  and  21 . The  drag  coefficient  behavior  clearly  demonstrates  the 
virtual  mass  effect  before  decaying  to  zero. 

4. 3 Non-Lifting  Flow  Past  an  Ellipse 

Inviscid  calculations  for  flow  past  a series  of  ellipses  of  various  thick- 
ness ratios  and  at  various  free  stream  Mach  numbers  have  been  carried  out  in  a 
manner  similar  to  that  for  the  cylinder  cases  described  above.  Sample  results 
are  presented  in  Figures  22  and  23  . 

The  family  of  ellipses  was  defined  by 

b = [l  + e sin2e]"1/2 

where  the  parameter  e is  related  to  the  body  thickness  ratio  h by 
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FIGURE  22.  ELLIPSE:  SURFACE  PRESSURE  DISTRIBUTION 


Catherall  - Sells  Solution 


Incompressible  Solution 
Time  Dependent  Solutior 


FIGURE  23.  ELLIPSE:  NEAR  FIELD  PRESSURE  DISTRIBUTION 
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The  transformation  function  x was  chosen  as 
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The  transformation  function  Y was  the  same  as  that  used  for  the  cylinder  flow 
field  computations  described  above  with  the  parameter  a determined  in  the  same 


Figure  22  compares  the  calculated  surface  pressure  distribution  for  a 
non- lifting  ellipse  of  thickness  ratio  h = 0.6  at  a free  stream  Mach  number  of 
n* 1 with  the  exact  incompressible  solution  and  the  Catherall-Sells  solution. 
Reasonable  accuracy  is  again  obtained  with  only  21  upper  surface  points  and 
30  points  in  the  radial  direction  even  though  the  body  shape  is  slightly  more 
complex  than  the  circular  cylinder.  The  angular  stretching  transformation  X 
defined  above  decreases  the  angular  spacing  in  the  vicinity  of  the  leading 
and  trailing  edges.  Figure  23  shows  the  calculated  pressure  distribution  near 
the  body  compared  with  the  incompressible  and  Catherall-Sells  solutions. 

4.4  Non-Lifting  Flow  Past  a Symmetric  Joukowsky  Airfoil 

Calculations  for  inviscid  non- lifting  flow  past  several  symmetric  Joukowsky 
airfoils  of  different  thickness  ratios  and  at  different  free  stream  Mach  numbers 
have  been  carried  out.  Sample  results  are  presented  in  Figures  24  through  32. 
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FIGURE  24.  JOUKOWSKY  AIRFOIL:  SURFACE  PRESSURE  DISTRIBUTION  (INVISCID) 


FIGURE  25.  JOUKOWSKY  AIRFOIL:  SURFACE  VELOCITY  DISTRIBUTION  (INVISCID) 


FiGURE  26.  JOUKOWSKY  AIRFOIL:  TIME  VARIATION  Of  STAGNATION, 
MINIMUM  AND  TRAILING  EDGE  PRESSURES  (INVISCID) 


FIGURE  28.  JOUKOWSKY  AIRFOIL:  SURFACE  VELOCITY  DISTRIBUTION  (INVISCID) 


FIGURE  30.  JOUKOWSKY  AIRFOIL:  SURFACE  PRESSURE  DISTRIBUTION  (INVISCID) 


FIGURE  31.  JOUKOWSKY  AIRFOIL:  SURFACE  VELOCITY  DISTRIBUTION  (INVISCID) 


FIGURE  32  JOUKOWSKY  AIRFOIL:  TIME  VARIATION  OF  STAGNATION. 

MINIMUM  AND  TRAILING  EDGE  PRESSURES  (INVISCID) 


The  transformation  functions  X and  Y chosen  for  this  family  of  body 
shapes  are  defined  by 
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The  quantities  p,  cf>  and  S(<(),t)  are  related  to  the  Joukowsky  transformation 
referred  to  earlier  in  Equations  (7)  through  (10). 

The  transformation  parameter  x,  an  input  quantity,  exerts  considerable 
influence  on  the  angular  distribution  of  body  grid  points  as  can  be  seen  from 
Figure  33.  Its  primary  effect  is  to  pull  points  from  the  trailing  edge  region 
and  cluster  them  near  the  surface  inflection  point  occurring  at  approximately 
mid-chord  while  leaving  the  nose  point  distribution  relatively  unchanged. 

The  arbitrary  functions  a.  and  6 were  chosen  as 

i(l>,t)  = t(t)  + 4 ( 1 1 s i n 

and 

o ( $,  t ) - FJt)  + '(t)sin  ’■} 

where  C,  ip , £ and  6 are  determined  to  control  the  radial  grid  spacing  at 
several  critical  locations.  Initially,  C and  £ are  determined  such  that  in 
the  physical  space  the  first  grid  interval  downstream  of  the  trailing  edge 
and  the  last  interval  just  inside  the  outer  wave  along  the  trailing  edge 
streamline  are  the  same  as  the  distance  between  the  trailing  edge  location  and 
the  nearest  adjacent  surface  grid  point.  This  grid  distance,  denoted  by  d ^ , 
is  typically  the  smallest  in  the  flow  field  and  is  used  in  the  CFL  rule  (Equation 
(103))  to  determine  the  allowable  time  step  AT.  Equal  upstream  and  downstream  soacine 
at  the  trailing  edge  is  advantageous  when  implementing  the  trailing  edge 
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analysis  of  Section  2.3.3  in  the  numerical  computations.  For  later  times, 
C and  £ are  assumed  to  vary  linearly  according  to 


C = + C 

wliere  the  proportionality  factor  Sf  is  an  arbitrary  input  quantity  (e.g., 

S = ~0.2)  and  C is  determined  from  the  initial  values  of  £ and  £.  When  the 
physical  space  begins  to  expand  (t  > t ) , the  upstream  and  downstream  trail- 
ing edge  grid  spacings  are  maintained  at  their  initial  values  and  this  condi- 
tion, together  with  the  above  linear  relationship,  is  sufficient  to  determine 
<;  and  £ at  all  later  times. 

The  grid  control  scheme  described  above  may  be  expressed  mathematically 
in  the  following  manner.  The  radial  coordinates  of  the  point  immediately 
downstream  of  the  trailing  edge  and  the  point  just  inside  the  outer  wave  on 
the  trailing  edge  streamline  are  initially 


r,  = b + d . 

2 min 


r » s(0,t  ) - d . 
n o mm 


where  d . is  determined  from  the  transformation  function  X.  From  the  Joukow- 
nun 

sky  transformation,  the  corresponding  values  of  p are 
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where  b denotes  the  trailing  edge  coordinate.  The  initial  values  of  ",  and  t, 
are  determined  to  satisfy 
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and 
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After  the  physical  space  begins  to  expand,  is  held  constant  and  £ and  £ 
are  determined  to  satisfy 


P2  = c + [s (0, t) 
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and 
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The  time  derivatives  of  C and  £ are  obtained  from  these  relations  as 


d£ 

dt 
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and 


d£  = d£ 
dt  C dt 

Similar  conditions  are  imposed  along  the  constant  X coordinate  line 

originating  from  the  first  surface  grid  point  upstream  of  the  body  surface 

inflection  point  where  a defined  in  Equations  (43)  vanishes.  That  is,  at 

time  zero  the  outward  spacing  at  the  body  surface  and  at  the  outer  wave  along 

this  grid  line  are  set  equal  to  d . to  determine  initial  values  of  ill  and  6. 

° min 

A linear  relationship  of  the  form 

tb  = S <$  + C. 

* ip  4) 

is  assumed,  where  S is  an  input  variable  (e.g.,  S = -.5)  and  C is  deter- 

4 ip  ip 

mined  from  the  initial  values  of  ip  and  6.  After  the  physical  space  begins  to 
expand,  the  outward  spacing  near  the  inflection  point  is  kept  equal  to  so 

that  ip  and  6 can  be  determined  at  all  later  times. 
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The  derivatives  of  the  transformation  functions  X and  V are  given  by 
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The  derivatives  of  p and  <j>  are  obtained  from  the  Joukowsky  transformation  as 
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Straightforward  differentiation  of  the  above  expressions  yields  the  derivatives 
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Figures  24  and  25  compare  the  calculated  surface  pressure  and  velocity 
distributions  for  an  airfoil  thickness  ratio  h = 0.118  at  a free  stream  Mach 
number  of  0.2  with  the  incompressible  solution  and  the  Catherall-Sells  solu- 
tion. Some  off-body  points  along  the  upstream  and  downstream  symmetry  stream- 
lines are  also  shown  to  indicate  the  extent  of  the  steady  flow  region  about 
the  airfoil.  Good  agreement  for  this  complex  shape  is  obtained  with  51  surface 
points  and  40  points  in  the  radial  direction.  The  accuracy  of  the  singular 
trailing  edge  analysis  of  Section  2. 3. 3.1  is  also  seen  to  be  quite  good.  The 
effectiveness  of  the  angular  transformation  function  X can  be  seen  from  the 
decreased  grid  spacing  in  the  neighborhood  of  the  surface  inflection  point. 

The  time  variation  of  stagnation,  minimum  and  trailing  edge  pressures  is  shown 
in  Figure  26  and  indicates  when  steady  state  conditions  around  the  body  are 
first  approached.  The  reason  for  this  particular  choice  of  thickness  ratio 
is  the  existence  of  experimental  data  (Reference  15),  including  boundary  layer 
measurements,  with  which  the  viscous  TDp  calculations  may  be  compared. 
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Figures  27,  28  and  29  show  comparisons  for  the  same  flow  conditions  and 
a much  larger  number  of  grid  points  (i.e.,  70  surface  points  and  50  radial 
points).  Increased  accuracy  is  obtained  especially  in  the  trailing  edpa 
region  and  the  surface  velocity  prediction.  Figures  30,  31  and  32  present 

I 

comparisons  for  an  airfoil  thickness  ratio  h = 0.10  and  a free  stream  Mach 
number  of  0.5.  These  predictions  were  made  with  70  surface  points  and  35  radial 
points.  The  accuracy  of  the  calculated  values  is  seen  to  be  reasonably  good. 


5.  VISCOUS  SAMPLE  CASES 


Viscous  calculations  were  carried  out  for  a non-lifting  Joukowsky  airfoil 
with  a thickness  ratio  of  0.118  to  permit  a comparison  with  the  test  results 
of  Reference  15.  These  test  results  include  detailed  surveys  of  the  viscous 
flow  in  the  boundary  layer  and  in  the  wake  in  the  neighborhood  of  the  trailing 
edge  suitable  for  evaluating  a viscous  calculation  method.  The  test  condi- 
tions were  essentially  incompressible  (M  = 0.036)  with  measurements  recorded 
for  both  a smooth  surface  airfoil  and  the  airfoil  with  a boundary  layer  transi- 
tion strip.  The  Reynolds  number  was  relatively  small  (Re  = 420,000)  having  the 
advantage  of  increasing  the  viscous  displacement  effects  but  adding  the  com- 
plication of  partly  turbulent  flow.  The  data  for  the  smooth  airfoil  were 
selected  for  comparison  purposes  because  the  TDC  program  has  no  provision  for 
representing  boundary  layer  thickening  due  to  a transition  strip. 

Inviscid  flow  and  boundary  layer  calculations  were  also  performed  for  com- 
parison with  the  test  data.  Figure  34  compares  the  experimental  pressure 
distribution  with  the  inviscid  flow  solutions  obtained  by  the  Catherall-Sells 
method  and  the  TDC  program  with  the  70x50  grid  discussed  previously  in  Section 
4.4.  The  inviscid  calculations  were  carried  out  for  M =0.2  because  of  the  more 

QO 

rapid  convergence  of  the  TDC  program  to  the  local  steady  state  condition  at 
the  higher  Mach  number.  The  compressibility  effects  resulting  from  the  dif- 
ference from  the  test  Mach  number  are  not  large  enough  to  affect  the  compari- 
son (see  Figures  27  and  28).  The  test  data  match  the  inviscid  calculations 
very  well  over  the  forward  portion  of  the  airfoil  (%  chord  <■  40)  , indicating 
good  quality  measurements  and  small  viscous  effects.  The  viscous  effects  are 
evident  as  a rounding  of  the  pressure  peak  at  the  trailing  edge  with  a trailing 
edge  Cp  about  halt  the  inviscid  value  and  with  Cp  values  in  the  wake  somewhat 
smaller  than  the  inviscid  values.  The  close  agreement  of  the  TDC  and  Catherall- 
oclls  calculations  near  the  trailing  edge  indicates  that  viscous  effects  of 
this  magnitude  should  be  predictable  by  the  time  dependent  method  of  computa- 
tion. 

Results  of  boundary  layer  surveys  aft  of  mid-chord  are  given  in  Reference 
15.  These  are  compared  with  results  of  the  Cebeci  finite-difference  boundary 
layer  method  (Reference  5)  in  Figures  35  through  38.  The  boundary  layer  cal- 
culation was  carried  out  with  the  inviscid  pressure  distribution  from  the 
Catherall-Sells  mpthod  and  also  with  the  measured  pressure  distribution  aft 
of  21.2%  chord.  Reference  15  also  gives  measurements  of  transition  location 
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FIGURE  34.  VISCOUS  EFFECT  ON  PRESSURE  DISTRIBUTION:  1 1 .8%  JOUKOWSKY  AIRFOIL 
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FIGURE  35.  COMPARISON  OF  BOUNDARY  LAYERS  ON  JOUKOWSKY  AIRFOIL: 
SHAPE  FACTOR 


FIGURE  36.  COMPARISON  OF  BOUNDARY  LAYERS  ON  JOUKOWSKY  AIRFOIL: 
MOMENTUM  THICKNESS 
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FIGURE  37.  COMPARISON  OF  BOUNDARY  LAYERS  ON  JOUKOWSKY  AIRFOIL: 
DISPLACEMENT  THICKNESS 
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FIGURE  38.  COMPARISON  OF  BOUNDARY  LAYERS  ON  JOUKOWSKY  AIRFOIL: 
BOUNDARY  LAYER  THICKNESS 
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for  the  smooth  airfoil  obtained  by  a smoke  filament  technique  which  indicate 
transition  near  mid-chord.  This  is  aft  of  the  laminar  separation  point  at  -8.0% 
chord  obtained  by  extrapolating  to  zero  the  laminar  skin  friction  results  from 
the  Cebeci  method.  At  such  a low  Reynolds  number  the  transition  process  should 
occur  gradually  over  a substantial  distance.  It  is  not  known  exactly  what  point 
within  the  transition  region  is  indicated  by  the  smoke  filament  technique, 
so  that  transition  may  have  started  ahead  of  the  measured  location.  Laminar 
separation  may  have  started  transition,  although  the  smoke  filament  did  not 
show  a separation  bubble.  Also,  the  separation  location  could  have  been 
changed  from  that  predicted  by  the  Cebeci  calculation  by  small  variations 
in  pressure  gradient  due  to  model  inaccuracies.  The  Cebeci  program  was  also 
used  to  estimate  a transition  location  by  repeating  the  boundary  layer  cal- 
culation for  several  assumed  transition  locations  and  selecting  that  transi- 
tion point  which  gave  the  best  agreement  with  the  measured  values  of  both 
shape  factor  and  momentum  thickness.  The  transition  point  location  is  an 
important  input  quantity  to  the  TDC  program  because  the  turbulence  equations 
employed  do  not  model  the  transition  process.  Abrupt  transition  is  assumed 
with  only  the  laminar  terms  in  the  mean-flow  equations  entering  the  calcula- 
tions upstream  of  some  preselected  boundary  in  the  computational  space, 
while  the  additional  turbulence  terms  and  turbulence  equations  are  employed 
beyond  this  boundary. 

The  Cebeci  method  has  options  for  either  gradual  or  abrupt  transition. 

It  was  felt  that  the  gradual  transition  option  would  better  represent  the 
actual  flow  conditions,  but  that  the  instantaneous  option  would  better  match 
the  results  from  a TDC  program  calculation.  The  best  match  with  the  test 
data  was  obtained  with  the  start  of  gradual  transition  at  -23.9%  chord  and 
with  the  point  of  instantaneous  transition  at  -9.6 % chord.  Figure  35  shows 
that  the  gradual  transition  calculation  agrees  better  with  the  measured 
shape  factor  data  indicating  that  transition  was  probably  not  complete  at 
the  first  two  experimental  measurement  stations.  However,  except  for  the 
first  station,  the  momentum  thickness  6 and  displacement  thickness  6*  from 
both  options  are  in  excellent  agreement  with  the  measurements.  This  implies 
that  either  abrupt  or  gradual  modelling  of  transition  will  give  an  equally 
good  prediction  of  the  viscous  interaction  effects  near  the  trailing  edge  if 
the  transition  location  is  adjusted  accordingly. 

The  rapid  drop  of  the  shape  factor  H a<.  the  point  of  instantaneous  transi- 
tion is  notable.  In  response  to  the  sudden  application  of  turbulent  eddy 
viscosity,  the  value  of  H at  this  point  has  completed  most  of  the  change  to 
an  equilibrium  turbulent  level  of  1.6  from  the  laminar  value  of  3.2  at  the 


86 


previous  station.  Combination  with  the  smooth  variation  of  momentum  thickness 
shown  in  Figure  36  results  in  a sharp  peak  in  the  displacement  thickness  at 
the  transition  point  (see  Figure  37).  This  might  be  expected  to  produce  larger 
viscous  effects  on  the  pressure  distribution  than  the  smoother  displacement 
shape  obtained  from  the  gradual  transition  calculation.  However,  this  dif- 
ference in  displacement  thickness  would  be  reduced  by  a fully  coupled  viscous- 
inviscid  solution. 

The  momentum  thickness  and  displacement  thickness  which  were  calculated 
from  the  inviscid  pressure  distribution  are  in  close  agreement  with  the  mea- 
sured values  near  the  trailing  edge.  However,  the  good  agreement  seems  to  be 
coincidental,  as  the  thicknesses  calculated  from  the  experimental  pressure  dis- 
tribution are  significantly  lower  than  the  measured  values.  The  experimental 
flow  field  data  show  a drop  in  pressure  outward  across  the  boundary  layer  at 
the  trailing  edge  so  that  the  disagreement  with  the  measured  thicknesses  would 
be  even  larger  if  the  pressure  at  the  boundary  layer  edge  were  used  in  the 
calculation.  It  seems  apparent  that  this  test  case  has  viscous-inviscid 
interaction  effects  near  the  trailing  edge  which  are  not  calculated  correctly 
by  the  boundary  layer  approximation,  but  should  be  calculated  correctly  by 
solving  the  turbulent  Navier-Stokes  equations  by  the  TDC  method. 

The  boundary  layer  calculations  were  used  as  a guide  in  selecting  the  grid 
geometry  shown  in  Figure  39  for  the  TDC  viscous  calculations.  A grid  with  70 
points  along  the  body  gave  good  accuracy  for  the  inviscid  calculations  (see 
Figure  34);  this  surface  point  distribution  was  therefore  used  for  the  viscous 
calculations.  A total  of  60  points  in  the  direction  normal  to  the  body  was 
used  to  provide  closer  spacing  in  the  viscous  region.  The  outer  boundary  for 
the  viscous  calculations  was  chosen  to  lie  well  outside  the  boundary  layer 
and  wake  so  that  the  viscous  terms  in  the  flow  equations  would  be  negligible 
at  the  boundary.  The  edge  of  the  viscous  flow  region  is  not  a well  defined 
boundary.  Several  estimates  of  the  boundary  layer  and  wake  thickness  S for 
the  selected  test  case  are  shown  in  Figure  38.  The  experimental  values  are 
based  on  a visual  inspection  of  the  measured  velocity  profiles,  and  probably 
correspond  to  a velocity  defect  of  one-half  or  one  percent  of  the  inviscid 
edge  velocity.  The  values  from  the  Cebeci  boundary  layer  calculations  are 
considerably  larger.  They  correspond  to  a numerical  test  which  selects  the 
outer  edge  of  the  boundary  layer  as  that  point  where  the  normal  velocity 
gradient  is  about  10  ^ of  the  surface  value.  The  selection  of  the  fifteenth 
grid  line  out  from  the  airfoil  (labeled  J = 15  in  Figure  39)  gave  a viscous 
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region  about  three  times  as  thick  as  the  Cebeci  limit.  The  grid  line  labeled 
I = 65  in  Figure  39  was  selected  as  the  outer  viscous  boundary  in  the  wake 
region. 

The  grid  line  labeled  I = 32  in  Figure  39  was  selected  as  the  upstream 
boundary  for  the  turbulent  calculations.  This  is  the  first  grid  line  upstream 
of  the  surface  point  corresponding  to  the  transition  location  determined  from 
the  Cebeci  instantaneous  transition  calculations.  Based  on  the  results  of  the 
turbulent  sample  case  presented  later,  this  may  not  have  been  the  right  choice. 

The  Cebeci  method  uses  an  empirical  mixing  length  formula  for  eddy  viscosity; 
the  large  increase  in  eddy  viscosity  at  the  instantaneous  transition  point 
causes  a rapid  change  in  the  boundary  layer  parameters,  as  was  shown.  Be- 
cause of  the  zero  upstream  boundary  condition  on  K and  E in  the  TDC  flow 
equations,  it  takes  some  distance  for  turbulence  to  be  produced  and  diffuse 
through  the  boundary  layer  to  develop  an  eddy  viscosity  corresponding  to  that 
used  in  the  Cebeci  method.  The  boundary  of  the  turbulent  region  could  possibly 
be  moved  a short  distance  upstream  to  allow  for  this.  In  the  test  case  that 
was  calculated,  this  eddy  viscosity  difference  is  important,  because  without 
fully  developed  turbulence  the  boundary  layer  may  separate  near  the  laminar 
separation  point  a short  distance  downstream  of  the  selected  turbulence 
boundary. 

.lore  testing  of  the  program  will  be  required  to  select  a method  for  prop- 
erly simulating  transition.  One  possibility  is  to  use  an  "equilibrium"  value 
of  K on  the  upstream  boundary  rather  than  K = 0,  although  some  analysis  would 
be  required  to  determine  what  value  to  use  in  a time  dependent  calculation. 

Another  possibility  would  be  to  multiply  the  eddy  viscosity  u by  an  inter- 
mittency  factor  to  simulate  the  gradual  transition  process  as  the  Cebeci  pro- 
gram does.  Neither  of  these  approaches  has  been  tested. 

The  "law-of-the-wall"  boundary  condition  for  turbulent  flow  developed  in 
Section  3.3  requires  that  the  first  grid  point  out  from  the  surface  be  in  the 
"law-of-the-wall"  region  of  the  boundary  layer,  where  the  "law-of-the-wall"  velocity 
profile  is  valid  and  the  turbulent  shear  is  approximately  equal  to  the  wall 
shear.  Transformation  parameters  were  selected  to  meet  this  requirement  near 
mid-chord  where  the  normal  spacing  of  the  grid  lines  is  controlled.  The  velocity 
profiles  from  the  Cebeci  boundary  layer  calculatior  with  instantaneous  transi- 
tion were  used  as  a guide  and  are  shown  in  Figure  40.  The  turbulent  profile 
at  mid-chord  shows  an  approximate  "law-of-the-wall"  region  between  y ■ 20 
and  v+  = 80.  The  selected  transformation  parameters  gave  the  grid  point 
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FIGURE  40.  JOUKOWSKY  AIRFOIL:  BOUNDARY  LAYER  VELOCITY  PROFILES  IN 
"LAW-OF-THE  WALL"  COORDINATES 


locations  shown.  The  first  grid  point  out  from  the  wall  (J  = 2)  is  at  a 

distance  y/c  = 0.002  while  the  value  of  y+  is  47,  well  within  the  "law-of- 

the-wall"  region.  Farther  back  on  the  airfoil,  the  upper  y+  limit  of  the 

"law-of-the-wall"  region  is  even  larger.  The  laminar  profile  at  -39%  chord 

is  located  near  the  pressure  peak  and  approximates  a flat  plate  or  Blasius 

profile  having  a nearly  linear  variation  (u  + = y+)  out  to  about  one-third 

of  the  boundary  layer  thickness.  The  transitional  profile  at  -9.6%  chord 

does  not  have  a "law-of-the-wall"  region  so  that  the  slip  velocity  u+  de- 

w 

rived  from  the  "law-of-the-wall"  boundary  condition  will  be  less  accurate 
in  the  short  distance  during  which  a "law-of-the-wall"  region  is  developing. 

The  assumption  that  the  shear  at  the  first  point  outward  from  the  sur- 
face equals  the  wall  shear  used  to  determine  the  wall  boundary  conditions 
given  by  Equations  122  is  examined  in  Figure  41.  The  shear  variation  at 
mid-chord  obtained  from  the  Cebeci  calculation  is  shown  and  the  approximation 
is  seen  to  be  valid  out  to  approximately  y+  = 80.  The  first  grid  point  (J  = 

2)  lies  well  within  this  region.  For  the  transitional  profile,  the  shear 
near  the  first  grid  point  is  larger  than  the  wall  value  due  to  the  larger 
velocity  gradients  shown  in  Figure  40.  This  will  tend  to  increase  the  wall 
values  of  K and  E predicted  by  Equations  122  in  the  transitional  region.  How- 
ever, the  lag  in  production  of  turbulent  energy  in  this  region  inherent  in 
the  TDC  calculations  will  have  an  opposite  effect. 

The  analysis  of  the  boundary  layer  calculations  indicates  that  the 
selected  grid  geometry  gives  an  inner  grid  point  within  the  "law-of-the-wall" 
region  at  mid-chord.  For  larger  Reynolds  numbers  or  locations  farther  down- 
stream, the  outer  boundary  of  the  "law-of-the-wall"  region  becomes  larger  in 
"law-of-the-wall"  coordinates  but  approaches  a nearly  constant  ratio  to  the 
momentum  thickness.  Figure  42  shows  an  estimate  of  the  outer  limit  yjj  of 
the  "law-of-the-wall"  region  from  a correlation  of  experimental  velocity 
profiles  with  shape  factor  H (Reference  16).  For  near  flat-plate  shape  factors 
(1.3  < H < 1.5)  the  outer  limit  is  approximately  the  same  as  the  momentum 
thickness.  For  boundary  layers  approaching  separation  (H  > 2.0),  the 
limit  is  reduced  to  yn/9  = 0.4.  This  indicates  that  a grid  spacing  near  the 
wall  of  y/c  = 0.002  at  mid-chord  should  give  an  acceptable  inner  grid  point 
location  at  larger  Reynolds  number  and  at  angles  of  attack  approaching  trailing 
edge  separation. 
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Shear  at  Mid-Chord 


5.1  Laminar  Flow 


The  laminar  flow  calculation  was  carried  out  to  a non-dimensional  time 
of  T = 4.5,  at  which  point  the  inviscid  calculation  has  essentially  reached 
steady  state  (see  Figure  29).  The  number  of  time  steps  required  to  reach  this 
time  point  was  about  50  percent  larger  than  for  the  inviscid  calculation  be- 
cause the  minimum  grid  spacing  d was  decreased  by  approximately  one  third. 

min 

The  allowable  time  step  size  given  in  Equation  103  is  directly  proportional 
to  the  minimum  grid  spacing. 

The  time  variation  of  stagnation,  minimum  and  trailing  edge  pressures 
for  laminar  flow  is  shown  in  Figure  43  and  compared  with  the  inviscid  results. 
The  variations  of  stagnation  and  minimum  pressure  are  very  similar  to  the  in- 
viscid case  as  expected,  as  the  laminar  boundary  layer  near  the  nose  of  the 
airfoil  is  very  thin.  There  is  a shift  in  pressure  levels  giving  an  error  in 
stagnation  pressure  several  times  larger  than  for  the  inviscid  case.  The  grid 
spacing  near  the  nose  is  about  the  same  as  for  the  inviscid  case  so  that  the 
difference  is  probably  the  result  of  the  zero  slip  velocity  boundary  condition 
combined  with  a grid  spacing  which  places  several  grid  points  in  the  first  row 
outward  from  the  body  outside  the  laminar  boundary  layer.  For  time  dependent 
calculations,  it  is  not  practical  to  use  a fine  grid  spacing  near  the  stagnation 
point  with  grid  points  in  the  constant  shear  layer  near  the  wall.  It  would  be 
desirable  to  develop  an  improved  numerical  scheme  which  would  permit  a rela- 
tively coarse  grid  while  giving  reasonably  accurate  results. 

The  minimum  pressure  reaches  a nearly  constant  level  near  T = 2 similar 
to  the  inviscid  result  and  is  in  fair  agreement  with  the  Catherall-Sells  steady 
state  value.  At  T = 4.5,  the  minimum  pressure  increases  to  a less  negative 
value  and  is  no  longer  steady  with  time.  The  cause  of  this  change  is  related 
to  the  surface  pressure  distributions  shown  in  Figure  44  and  the  velocity 
profiles  shown  in  Figure  45.  The  pressure  is  less  negative  over  the  forward 
part  of  the  airfoil  than  the  steady  state  inviscid  solution  and  much  more 
negative  over  the  rear  part.  The  velocity  profiles  of  Figure  45  show  a 
change  from  a relatively  thick  boundary  layer  type  profile  at  T = 1.9  to  a 
very  thick  shear  layer  type  profile  at  T = 4.5.  The  boundary  layer  type  pro- 
files (over  the  forward  part  of  the  airfoil  surface)  increase  gradually  in 
thickness  in  the  early  part  of  the  time  history.  The  thickening  of  the 
boundary  layer  profiles  may  represent  the  initial  stages  of  laminar  separa- 
tion as  predicted  by  the  boundary  layer  program,  but  the  final  shear  layer 
profile  Is  meaningless  since  the  computation  near  the  trailing  edge  is  going 
unstable. 


94 


FIGURE  43.  JOUKOWSKY  AIRFOIL:  TIME  VARIATION  OF  STAGNATION, 
MINIMUM  AND  TRAILING  EDGE  PRESSURES  (VISCOUS) 


FIGURE  44.  JOUKOWSKY  AIRFOIL:  SURFACE  PRESSURE  DISTRIBUTION  (LAMINAR) 


The  time  variation  of  trailing  edge  pressure  shown  in  Figure  43  is  much 
different  from  the  inviscid  case.  It  is  nearly  constant  and  approximately 
zero  at  the  start  of  the  motion  followed  by  a drop  to  a constant  negative 
value.  It  should  be  expected  to  follow  the  inviscid  curve  more  closely, 
especially  at  early  times  when  the  boundary  layer  is  thinner.  The  differ- 
ence seems  to  be  due  to  trailing  edge  numerical  problems  in  the  evaluation 

of  jlP  . An  additional  indication  of  this  is  the  computational  noise  near 
3T 

the  trailing  edge  shown  in  Figure  44  which  grows  in  amplitude  with  distance 
as  the  trailing  edge  is  approached.  Similar  computational  noise  was  expe- 
rienced in  the  inviscid  calculations  before  more  accurate  numerical  tech- 
niques to  evaluate  the  trailing  edge  pressure  were  developed.  When  the  im- 
proved techniques  were  added  to  the  program,  the  computational  noise  de- 
creased dramatically  giving  the  smooth  inviscid  results  shown  in  Section  4.4. 

A local  mesh  refinement  may  also  improve  the  stability  near  the  trailing  edge. 
5. 2 Turbulent  Flow 

The  turbulent  calculations  were  carried  out  to  a non-dimensional  time  of 
only  T = 2 because  of  computational  cost.  At  this  time  the  rapid  transients 
similar  to  those  of  the  inviscid  calculation  have  subsided  and  only  gradual 
changes  occur  in  the  subsequent  approach  to  steady  state.  The  turbulent  cal- 
culations behave  in  a manner  very  similar  to  the  laminar  results  in  many 
respects,  and  the  later  time-wise  flow  development  should  show  some  of  the 
same  trends. 

The  time  variation  of  stagnation,  minimum  and  trailing  edge  pressures 
for  turbulent  flow  is  shown  in  Figure  43  and  compared  to  the  inviscid  and 
laminar  results.  The  stagnation  and  minimum  pressure  variations  are  almost 
identical  to  the  laminar  case  as  expected,  since  the  forward  part  of  the 
viscous  flow  is  laminar.  The  minimum  pressure  variation  is  more  irregular 
for  the  laminar  case,  because  the  turbulent  case  has  reduced  computational 
noise  over  the  aft  portion  of  the  airfoil  as  shown  in  Figure  46.  This  is 
probably  due  to  the  damping  effect  of  turbulence  which  acts  to  reduce  the 
noise  magnitude  considerably. 

In  Figure  45  the  velocity  profile  at  mid-chord  for  the  turbulent  case  at 
T = 2 is  compared  to  the  laminar  result  at  the  closest  time  for  which  data 
were  available  (i.e.,  T = 1.9).  The  profiles  are  very  similar  and,  allowing 
for  the  small  time  difference,  there  is  no  apparent  effect  of  turbulence  except 
for  the  slip  velocity  at  the  surface  in  the  turbulent  case.  Turbulence  might 
be  expected  to  have  a stronger  effect  on  the  velocity  profile,  considering  the 
rapid  change  in  displacement  thickness  6*  and  shape  factor  H aft  of  the 
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FIGURE  46.  JOUKOWSKY  AIRFOIL:  SURFACE  PRESSURE  DISTRIBUTION  (TURBULENT) 


transition  point  as  determined  by  the  Cebeci  method.  However,  it  was  noted 
before  that  the  gradual  production  of  turbulence  by  the  TDC  flow  equations 
downstream  of  the  boundary  of  turbulent  flow  does  not  correspond  to  the 


sudden  increase  in  eddy  viscosity  inherent  in  the  Cebeci  calculation.  The 
boundary  layer  at  mid-chord  may  be  in  the  early  stages  of  laminar  separa- 
tion in  both  the  laminar  and  turbulent  cases.  It  is  recommended  that  addi- 
tional turbulent  calculations  be  performed  with  the  boundary  of  turbulent 
flow  located  farther  upstream  to  avoid  the  tendency  toward  laminar  separa- 
tion. Any  improvements  in  the  laminar  boundary  layer  numerical  techniques 
would  likewise  apply  to  the  turbulent  case,  as  the  same  equations  are  used 
for  the  laminar  boundary  layer  portion  of  the  flow. 


6.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  inviscid  version  of  the  Time  Dependent  Computation  (TDC)  program  has 
been  used  to  calculate  the  flow  field  about  various  body  shapes.  For  flows 
with  lift,  introduction  of  circulation  into  the  far  field  has  been  found  to 
be  the  physically  correct  means  of  generating  lift  in  an  inviscid  flow.  The 
validity  of  the  method  of  characteristics  approach  to  correctly  predict  bound- 
ary conditions  at  a solid  surface  or  moving  discontinuity  has  also  been  demon- 
strated. A technique  for  handling  singular  boundary  points  such  as  cusped 
trailing  edges  has  been  presented  and  its  accuracy  verified  for  the  inviscid 
case  of  a symmetric  Joukowsky  airfoil. 

A two-equation  turbulence  model  was  selected  to  extend  the  TDC  program 
formulation  to  turbulent  flow.  The  additional  cost  of  viscous  calculations 
was  reduced  by  limiting  the  evaluation  of  viscous  terms  to  a region  encom- 
passing the  body  and  the  wake.  Additional  economy  was  realized  for  turbulent 
flow  cases  by  using  a "law-of-the-wall"  boundary  condition  in  the  turbulent 
region  near  the  surface  to  avoid  the  small  grid  spacing  required  to  define  the 
steep  gradients  existing  in  the  viscous  sublayer  region.  Analysis  of  boundary 
layer  calculations  showed  that  the  minimum  grid  spacing  required  by  the 
"law-of-the-wail"  boundary  condition  need  not  be  much  smaller  than  that  required 
for  an  accurate  inviscid  calculation. 

Viscous  flow  calculations  for  both  completely  laminar  and  partially  tur- 
bulent boundary  layer  and  wake  flow  were  carried  out  with  the  TDC  program  for 
a non-lifting,  symmetrical  Joukowsky  airfoil  at  a low  Reynolds  number.  The 
calculated  pressures  near  the  airfoil  nose  for  the  laminar  flow  case  were  not 
as  accurate  as  for  the  inviscid  case,  even  though  the  grid  point  spacing  was 
slightly  denser,  indicating  a need  for  improvement  in  the  numerical  techniques 
used  to  evaluate  the  laminar  viscous  boundary  condition.  Boundary  layer 
type  velocity  profiles  increased  gradually  in  thickness  in  the  early  part 
of  the  time  history.  At  later  times,  the  flow  field  results  developed  an 
instability  near  the  trailing  edge,  with  a large  amount  of  computational 
noise  resulting  from  the  inadequate  numerical  techniques  used  at  the  singular 
trailing  edge  point  for  laminar  flow. 

The  turbulent  flow  calculation  was  carried  out  with  the  turbulence 
limited  to  the  region  downstream  of  a "transition"  boundary.  This  boundary 
was  selected  as  the  transition  point  for  which  independent  boundary  layer 
calculations  produced  a boundary  layer  development  agreeing  with  the  exist- 
ing test  data  of  Reference  15.  The  addition  of  turbulence  effects  to  the 
calculation  did  not  significantly  alter  the  thickening  of  the  viscous  profile 
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as  observed  in  the  laminar  case.  The  transition  boundary  may  have  been 
located  too  short  a distance  upstream  of  the  predicted  laminar  separation 
point  to  produce  the  expected  turbulent  boundary  layer  development.  The 
turbulence  considerably  reduced  the  magnitude  of  the  computational  noise  in 
the  trailing  edge  region  although  not  sufficiently  to  give  accurate  results. 
This  indicates  the  trailing  edge  numerical  techniques  also  need  improvement 
for  turbulent  flow. 

It  is  recommended  that  improved  numerical  techniques  for  the  viscous 
boundary  conditions  be  developed  and  tested,  with  special  attention  given  to 
the  trailing  edge  singularity.  A better  method  of  simulating  transition 
should  also  be  developed.  These  improvements  are  needed  before  the  capa- 
bility of  the  TDC  method  to  predict  viscous/inviscid  interaction  effects  can 
be  adequately  determined. 
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